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Abstract. 

Numerical integration of the differential equations of light propagation in the 
Schwarzschild metric shows that in some situations relevant for practical observations 
the well-known post-Newtonian solution for light propagation has an error up to 16 
//as. The aim of this work is to demonstrate this fact, identify the reason for this 
error and to derive an analytical formula accurate at the level of 1 /ias as needed for 
high-accuracy astrometric projects (e.g., Gaia). 

An analytical post-post-Newtonian solution for the light propagation for both 
Cauchy and boundary problems is given for the Schwarzschild metric augmented by 
the PPN and post-linear parameters /?, 7 and e. Using analytical upper estimates 
of each term we investigate which post-post-Newtonian terms may play a role for an 
observer in the solar system at the level of 1 /xas and conclude that only one post- 
post-Newtonian term remains important for this numerical accuracy. In this way, an 
analytical solution for the boundary problem for light propagation is derived. That 
solution contains terms of both post-Newtonian and post-post-Newtonian order, but 
is valid for the given numerical level of 1 //as. The derived analytical solution has 
been verified using the results of a high-accuracy numerical integration of differential 
equations of light propagation and found to be correct at the level well below 1 /ias 
for arbitrary observer situated within the solar system. Furthermore, the origin of the 
post-post-Newtonian terms relevant for the microarcsecond accuracy is elucidated. We 
demonstrate that these terms result from an inadequate choice of the impact parameter 
in the standard post-Newtonian formulas. Introducing another impact parameter, that 
can be called "coordinate-independent" , we demonstrate that all these terms disappear 
from the formulas. 



PACS numbers: 95.10.Jk, 95.10.Ce, 95.30. Sf, 04.25.Nx, 04.80.Cc 



1. Introduction 

It is well known that adequate relativistic modeling is indispensable for the success 
of microarcsecond space astrometry. One of the most important relativistic effects 
for astrometric observations in the solar system is the gravitational light deflection. 
The largest contribution in the light deflection comes from the spherically symmetric 
(Schwarzschild) parts of the gravitational fields of each solar system body. Although 
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the planned astrometric satellites Gaia, SIM, etc. will not observe very close to the 
Sun, they can observe very close to the giant planets also producing significant light 
deflection. This poses the problem of modeling this light defiection with a numerical 
accuracy of better than 1 yuas. 

The exact differential equations of motion for a light ray in the Schwarzschild field 
can be solved numerically as well as analytically. However, the exact analytical solution 
is given in terms of elliptic integrals, implying numerical efforts comparable with direct 
numerical integration, so that approximate analytical solutions are usually used. In fact, 
the standard parametrized post-Newtonian (PPN) solution is sufficient in many cases 
and has been widely applied. So far, there was no doubt that the post-Newtonian order 
of approximation is sufficient for astrometric missions even up to microarcsecond level of 
accuracy, besides astrometric observations close to the edge of the Sun. However, a direct 
comparison reveals a deviation between the standard post-Newtonian approach and the 
high-accuracy numerical solution of the geodetic equations. In particular, we have found 
a difference of up to 16 fias in fight deflection for solar system objects observed close to 
giant planets. This error has triggered detailed numerical and analytical investigation 
of the problem. 

Usually, in the framework of general relativity or the PPN formalism analytical 
orders of smallness of various terms are considered. Here the role of small parameter 
is played by c~^, where c is the light velocity. Standard post-Newtonian and post- 
post-Newtonian solutions are derived by retaining terms of relevant analytical orders of 
magnitude. On the other hand, for practical calculations only numerical magnitudes of 
various terms are relevant. In this work we attempt to close this gap and combine the 
analytical parametrized post-post-Newtonian solution with exact analytical estimates of 
the numerical magnitudes of various terms. In this way we derive a compact analytical 
solution for light propagation where all terms are indeed relevant at the level of 1 //as. 
The derived analytical solution is then verified using high- accuracy numerical integration 
of the differential equations of light propagation and found to be correct at the level 
well below 1 /las. 

We use fairly standard notations: 

• G is the Newtonian constant of gravitation. 

• c is the velocity of light. 

• /3 and 7 are the parameters of the Parametrized Post-Newtonian (PPN) formalism 
which characterize possible deviation of the physical reality from general relativity 
theory = 7 = 1 in general relativity) . 

• Lower case Latin indices i, j, . . .take values 1, 2, 3. 

• Lower case Greek indices ji, v, .. .take values 0, 1, 2, 3. 

• Repeated indices imply the Einstein's summation irrespective of their positions (e.g. 
a' h' ^aH^ + + and a" b"" = 0° 6° + b^ + aH^ + b^). 
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• A dot over any quantity designates the total derivative with respect to the 

coordinate time of the corresponding reference system: e.g. a = —. 

at 

• The 3-dimensional coordinate quantities ( "3- vectors" ) referred to the spatial axes 
of the corresponding reference system are set in boldface: a = a\ 

• The absolute value (Euclidean norm) of a "3- vector" a is denoted as |a| or, simply, 
a and can be computed as a = \a\ = {a^ a^ + a'^ a'^ + a^ a^Y^'^. 

• The scalar product of any two "3-vectors" a and b with respect to the Euclidean 
metric 6ij is denoted by a ■ 6 and can be computed as a ■ 6 = 6ij a* IP = a^lf. 

• The vector product of any two "3-vectors" a and h is designated by a x 6 and can 
be computed as (a x 6)* = Sijk a^ 6^, where Sijk = {i — — k){k — i)/2 is the fully 
antisymmetric Levi-Civita symbol. 

• For any two vectors a and b, the angle between them is designated as 5{a,b). 
Clearly, for an angle between two vectors one has < 6{a, b) < tt. Angle 6{a, b) 

can be computed in many ways, for example, as 6(a, b) = arccos — — . 

a b 

This paper is a concise exposition of the work performed in the framework of the 
ESA project Gaia and published in a series of preprints [H] EH [121 ES]. The paper 
is organized as follows. In Section [2] we present the exact differential equations for 
the light propagation in the Schwarzschild field in harmonic gauge. High-accuracy 
numerical integrations of these equations are discussed in Section [31 In Section [H we 
discuss the standard post-Newtonian approximation and demonstrate its errors by direct 
comparison with numerical results. In Section [5] the analytical post-post-Newtonian 
solution for the light propagation is given. Section [61 is devoted to the boundary problem 
for the light propagation in post-post-Newtonian approximation. Investigations of the 
post-post-Newtonian terms in the formulas for the light deflection reveal that these 
terms can be divided into two groups: "regular" (those which can be estimated as 

2 

const - where m is the Schwarzschild radius of the deflecting body and d is the impact 
parameter) and "enhanced" (those which cannot be estimated like this and may become 
substantially larger than the "regular" terms). In Section [7] we clarify the physical origin 
of the "enhanced" post-post-Newtonian terms. The results are summarized in Section [HI 



2. Schwarzschild metric and null geodesies in harmonic coordinates 

We need a tool to calculate the real numerical accuracy of some analytical formulas 
for the light propagation in various situations. To this end, we consider the exact 
Schwarzschild metric and its null geodesies in harmonic gauge. Those exact differential 
equations for the null geodesies will be solved numerically with high accuracy (see below) 
and that numerical solution provides the required reference. 
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2.1. Metric tensor 
In harmonic gauge 



the components of the covariant metric tensor of the Schwarzschild solution are given 
by 

1 — a 
1 + a 

fi'oi = 0, 

g,^ = (l + af S,, + ^l^x'x\ (2) 



where 



a = — , (3) 

and the Schwarzschild radius of a body with mass M. The contravariant 

components of the metric read 
00 ^ l + a 
1 — a 

• ■ 1 1 ■ . 

(1 + a) 2; (1 + a) 

Considering that the determinant of the metric can be computed as 

g = -{l + a)\ (5) 

one can easily check that this metric satisfies the harmonic conditions ([1]). 

2.2. Christojfel symbols 

The Christoffel symbols of second kind are defined as 

"'^ 2 ^ \dx^ dx'- dx- J ■ 

Using ([2]) and (jlj) one gets 

_ ^ , 
°' x2 l-a2 ^ ' 
_ g l-a i 
(l + a)3'' ' 

r;. = ^ - ;^ (-^' + - ^ ^ -(7) 

All other Christoffel symbols vanish. 
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2.3. Isotropic condition 

The conditions that a photon follows an isotropic geodesic can be formulated as an 
equation for the four components of the coordinate velocity x°' of that photon: 

dx°' dx^ . , 

A being the canonical parameter, or 
2 1 

fi-oo + -9oii' + gij x'x^ =0, (9) 

where = dx'^/dt is the coordinate velocity of the photon. Eq. (j9]) is a first integral 
of motion for the differential equation for light propagation and must be valid for any 
point of an isotropic geodesic. Substituting the ansatz x = csfj,, where /i is a unit 
coordinate direction of light propagation (/x ■ /x = 1) and s = \x\/c, into ([9]) one gets for 
metric ([2]): 

s=\^ (l-a^ + -,{x.^y] (10) 
1 + a Y / 

This formula allows one to compute the absolute value of coordinate velocity of light in 
the chosen reference system if the position of the photon and the coordinate direction 
of its propagation fx are given. 



2.4^. Equation of isotropic geodesies 

The geodetic equations 

d'^x^ ^ dx" dx^ 
H)^ ^ ''^'dXliX^ 
can be re-parametrized by coordinate time t to give 



x' 



-c' T\,-2c x^ - V), x^ x^ + x^ [c V\, + 2 1% x^ + 1 r°, x^ i'^) 



(12) 



Substituting the Christoffel symbols one gets the differential equations for the light 
propagation in metric ([2]): 



1 — a 



i ■ i + a 



a f X ■ X 



1 -a2 



a 2 — a , ... , , 

x + 2—- -{x-x)x. (13) 

1 — a"^ 



Eq. ( ITU]) for the isotropic condition together with x ■ x = s"^ could be used to avoid 
the term containing x ■ x, but this does not simplify the equations and we prefer not to 
do this here. 



3. Numerical integration of the equations of light propagation 

Our goal is to integrate (IT^ numerically to get a solution for the trajectory of a light ray 
with an accuracy much higher than the goal accuracy of 1 /ias ~ 4.8 x 10~^^ radians. 
For these numerical integrations a simple FORTRAN 95 code using quadruple (128 
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bit) arithmetic has been written. Numerical integrator ODEX [7] has been adapted to 
the quadruple precision. ODEX is an extrapolation algorithm based on the explicit 
midpoint rule. It has automatic order selection, local accuracy control and dense 
output. Using forth and back integration to estimate the accuracy, each numerical 
integration is automatically checked to achieve a numerical accuracy of at least 10~^^ 
in the components of both position and velocity of the photon at each moment of time. 

The numerical integration is first used to solve the initial-value (Cauchy) problem 
for differential equations ( |T3l) . Eq. ( |TOl) should be used to choose the initial conditions. 
The problem of light propagation has thus only 5 degrees of freedom: 3 degrees of 
freedom correspond to the position of the photon and two other degrees of freedom 
correspond to the unit direction of light propagation (of course, in the Schwarzschild 
field with its symmetry one has also further integrals of motion, but here we ignore 
this; see Section [53] below). The absolute value of the coordinate hght velocity can 
be computed from (fTOj) . Fixing initial position of the photon x(to) and initial (unit) 
direction of propagation /j, one gets the initial velocity of the photon as function of fi 
and s computed for given /j, and x as given by (|TOl) : 

x{to) = Xo , 

x{to) = OS fi. (14) 

The numerical integration yields the position x and velocity x of the photon as function 
of time t. The dense output of ODEX allows one to obtain the position and velocity of 
the photon on a selected grid of moments of time. Eq. (fTOj) must hold for any moment 
of time as soon as it is satisfied by the initial conditions. Therefore, Eq. (ITOi) can also 
be used to check the accuracy of numerical integration. 

For the purposes of this work we need to have an accurate solution of two-value 
boundary problem. That is, a solution of Eq. ( fT3l) with boundary conditions 

x(to) = Xo, 

x{ti)=xi, (15) 

where Xq and Xi are two given constants, to is assumed to be fixed and ti is unknown 
and should be determined by solving (fT3|) . Instead of using some numerical methods 
to solve this boundary problem directly, we generate solutions of a family of boundary 
problems from our solution of the initial value problem (|T^ . Each intermediate result 
computed during the numerical integration with initial conditions ( fT4l) gives us a high- 
accuracy solution of the corresponding two- value boundary problem (fTSll : ti and Xi are 
simply taken from the numerical integration. 

As discussed in [10], the light propagation is characterized by three unit vectors 
(see Figured]): the coordinate direction n of light propagation at the point of reception 

the coordinate direction cr of light propagation for time going to minus infinity, 

cr = lim -x{t) , (17) 
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and the coordinate direction k from the point of hght emission to the point of reception 
R 

fc = — , R = Xi-xo. (18) 

In the following discussion we will compare predictions of various analytical models 
for n in the framework of the boundary problem ( |T5i) . The reference value for these 
comparisons can be computed using ( [T6|) and x(ti) from the numerical integration. The 
accuracy of this n computed from our numerical integrations is guaranteed to be of the 
order of 10~^^ radians and can be considered as exact for our purposes. 



4. The deficiency of the standard post-Newtonian approach 

Let us now demonstrate that the standard post-Newtonian formulas for the light 
propagation have too large numerical errors when compared to the accurate numerical 
solution of the geodetic equations described in the previous Section. 



4.1. Equations of the post- Newtonian approach 

The well-known equations of light propagation in first post-Newtonian approximation 
with PPN parameters have been discussed by many authors (see, for example, 
[201 [31 [iTj). Let us here summarize the standard post-Newtonian formulas. The 
differential equations for the light rays read (see also Section 15.1.41 below) 



X 



- (c' + jx'x') ^ + 2(l + 7)^^^^ + 0(c-2). (19) 

The analytical solution of ( !T9|) can be written in the form 

x{t) = a^pN + 0(c-^) , (20) 
a^pN = Xf) + c{t - to) a- + Ax{t) , (21) 

where 

/ / \ \ \ X -\- (T ■ X 

Ax{t) = — (1 + 7)m cr x (xq x cr) + cr log ■ 



(T ■ X Xq — (T ■ Xq/ Xo + (T ■ Xo 



(22) 
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Solution fl20l) - fl22l) satisfies the following initial conditions: 

x{to) = Xq , 

lim x{t) =c(T. (23) 

t— oo 

From (!20l) - (!22l) it is easy to derive the following expression for the unit tangent vector 
at the observer's position Xi for the boundary problem (|T5l) (the standard technique 
to do this is given, e.g. in |3] and used below in Section |6] in the post-post-Newtonian 
approximation) : 

np. = fc-(l + 7)m- , (24) 

where R and k are defined by ( |T8|) . and d = k x [xq x k) = k x [xi x k) is the impact 
parameter of the straight line connecting Xq and Xi. 

4-2. Comparison of the post- Newtonian formula and the numerical solution 

In order to investigate the accuracy of the standard post-Newtonian formulas we have 
compared the post-Newtonian predictions of the light deflection with the results of the 
numerical solution of geodetic equations. Here, we calculate the angle between the unit 
tangent vector ripN defined by (l2ll) and the vector n computed using (fT6l) from the 
numerical integration of (fT3!l . 

Having performed extensive tests, we have found that, in the real solar system, the 
error of for observations made by an observer situated in the vicinity of the Earth 
attains 16 /zas. These results are illustrated by Table [H and Figure O Table [U contains 
the parameters we have used in our numerical simulations as well as the maximal angular 
deviation between ripN and n in each set of simulations. We have performed simulations 
with different bodies of the solar systems, assuming that the minimal impact distance 
d is equal to the radius of the corresponding body, and the maximal distance Xi = \xi\ 
between the gravitating body and the observer is given by the maximal distance between 
the gravitating body and the Earth. The simulation shows that the error of UpN is 
generally increasing for larger xi and decreasing for larger d. The dependence of the 
error of rip^ for fixed d and Xi and increasing distance Xq between the gravitating 
body and the source is given on Figure [2] for the case of Jupiter, where minimal d and 
maximal xi (according to Table [1]) were used. Moreover, the error of npN is found to be 
proportional to which leads us to the necessity to deal with the post-post-Newtonian 
approximation for the light propagation. 



5. Analytical post-post-Newtonian solution 

The goal of this Section is to derive a rigorous analytical post-post-Newtonian solution 
for light propagation in the gravitational field of one spherically symmetric body in the 
framework of the PPN formalism extended by a non-linear parameter for the terms of 
order c~'^ in gij. The geodetic equation for the light ray in Schwarzschild metric can 
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Figure 2. The angle between ripN and n for Jupiter. The vector ripN is evahiated 
by means of the standard post-Newtonian formula ([M)) . while n is taken from the 
numerical integration as described in Section |31 Impact parameter d is taken to be the 
radius of Jupiter and the distance xi between Jupiter and the observer is 6 au. 
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Table 1. Numerical parameters of the Sun and giant planets are taken from 15], 
dmin is the minimal value of the impact parameter d that was used in the simulations. 
For each body, dmin is equal the radius of the body. For the Sun at 45° the impact 
parameter is computed as d = sin 45° x 1 au. x'^^^ is the maximal absolute value of 
the distance xi between the gravitating body and the observer that was used in the 
simulations. 5max is the maximal angle between ripN and n found in our numerical 
simulations. 





Sun 


Sun at 45° 


Jupiter 


Saturn 


Uranus 


Neptune 


TO = GM/c^ [m] 
d,„i„ [10^ m] 
x'l'^'' [au] 


1476.6 

696.0 

1 


1476.6 

105781.7 

1 


1.40987 
71.492 

6 


0.42215 

60.268 

11 


0.064473 

25.559 

21 


0.076067 

24.764 

31 




3187.8 


6.32 X lO-'' 


16.13 


4.42 


2.58 


5.84 



in principle be integrated exactly ^. However, such an analytical solution is given 
in terms of elliptic integrals and is not very suitable for massive calculations. Besides 
that, only the trajectory of the photon is readily available from the literature, but 
not the position and velocity of a photon as functions of time. Fortunately, in many 
cases of interest approximate solutions are sufficient. The standard way to solve the 
geodetic equation is the well-known post-Newtonian approximation scheme. Normally, 
in practical applications of relativistic light propagation, the first post-Newtonian 
solution is used. Post-post-Newtonian effects have been also sometimes considered 
[9l [13], but in a way which cannot be called self-consistent since no rigorous solution 
in the post-post-Newtonian approximation has been used. Such a rigorous post-post- 
Newtonian analytical solution for light propagation in the Schwarzschild metric has been 
derived in [21 [3] in general relativity in a class of gauges. However, the parametrization 
in [21 [3] does not allow one to consider alternative theories of gravity and therefore, 
a post-post-Newtonian solution for light propagation within the PPN formalism and 
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its extension to the second post-Newtonian approximation is not known. However, it 
is clearly necessary to have such a solution. Therefore, our goal is to generalize the 
post-post-Newtonian solution of [2] and to extend it for the boundary problem for light 
propagation. 

5.1. Differential equations of light propagation and their integral 

The first part of the problem is to derive the differential equations of light propagation 
with PPN and post-linear parameters. 

5.1.1. Metric tensor in the parametrized post-post-Newtonian approximation. 
Expanding metric ([2]) in powers of c~^, retaining only the terms relevant for the post- 
post-Newtonian solution for the light propagation, and introducing the PPN parameters 
(3 and 7 |2Qj and the post-linear parameter e one gets 

^00= -l + 2a-2Pa^ + 0{c~''), 
goi = 0, 

gij = 6ij + 2-fa6^J + e (^6ij + + 0{c-^) , (25) 

a being again defined by ([3]). In general relativity one has /3 = 7 = e = 1. Parameter e 
should be considered as a formal way to trace, in the following calculations, the terms 
coming from the terms in gij. No physical meaning of e is claimed here. However, 
this parameter is equivalent to parameter A of [HI |15[ [16] and parameter e of [6] . 

The corresponding contravariant components of metric tensor can be deduced from 
(125|1 and are given by 

^00= _ 1 _2a + 2(/3-2)a^ + C»(c-^), 

g^^ = 5., - 2 7 a 6,, + (^{4^' - e) 6,, - e a' + 0{c-'). (26) 

The determinant of metric tensor reads 

g = -l-2(37-l)a-2(/3 + 2e + 67(7-l))a2 + C>(c-''), (27) 

^ = l + (37-l)a+(2/3 + 4e-l + 37(7-2))a2 + 0(c-^). (28) 

Metric (125!) is obviously harmonic for 7 = /3 = e = 1 since the harmonic conditions ([1]) 
take the form 

9(^=^^°"^ 



^(vgs^ = (1 - ^) + ((1 + - 2/? - 2.) + O(c-O). (29) 
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5.1.2. Christoffel symbols. The Christoffel symbols of second kind defined by (E]) can 
be derived from metric f l25|l - fl26|) : 

r°oo = 0, (30) 

T\ = ^ + {l-^)^-^ + 0{c-% (31) 

r°.fc = 0, (32) 

^\o = ^-iP + l)^-^ + 0{c-% (33) 

= 0, (34) 
r'fcz = 7 {x' hi - 6ii - 6ik) 



i k l\ 2 

+ (2(e-7')x^5H-(e-272) (^x' 5u + 5,k) - 2 e^^^j ^ + 0(c-6). (35) 

5.1.3. Isotropic condition for the null geodetic. From now on, denote the coordinates 
of a photon, denote the spatial coordinates of the photon, and x = \x\ is the distance 
of the photon from the gravitating body that is situated at the origin of the used 
reference system. As it was discussed in Section 12.31 Eq. ^ allows one to compute the 
absolute value of coordinate velocity of light if the position of the photon x* and the 
unit coordinate direction of its propagation /i* (/x ■ /x = 1) are given. Using (125|1 for 
s = \x\/c one gets 

s = 1 - (1 + 7) a + i 1^-1 + 2/3 - e + 7 (2 + 37) - e (^^)'^ + C(c"^) • (36) 

5.1.4. Differential equations of light propagation. Inserting the Christoffel symbols 
fl5D]) - fl5Sl) into (fT^ . one gets the following equations of light propagation in post-post- 
Newtonian approximation 

ax , ^ axix ■ x\ 



X 



(c2 + 7i-i) — + 2(1 + 7) 



■ ■ ■ x^ 



+ 2 + 7) + (72 - e) (i . i)) ^ + 26 



a^ X „ a^ xix • xY 



+ 2 (2(1 - /3) + e - 2 7') -^^r + 0{c-') . (37) 



x^ 

Here, for estimating the analytical order of smallness of the terms we take into account 
that \x\ = 0{c). Using ( 136|) and x ■ x = c^ s^ one can simplify (1371) to get 

, . oO-^ r. . ax(x ■ x) 
x= - 1 + 7 c^ — + 2(1 + 7) — - 

+ 2c^ (/3 -6 + 27(1 + 7)) ^ + 2e 

+ 2 (2(1 - /3) + e - 2 7') ^^^^^ + 0(c-^) . (38) 
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5.1.5. Equations of light propagation with additional trace parameter a. For our 
purposes it is advantageous to have one more parameter that can be used to trace 
terms in the following calculations which come from the post-post-Newtonian terms in 
the equations of motion of a photon. We denote this parameter a and introduce it 
in the above equation simply as a factor for all the post-post-Newtonian terms in the 
right-hand side: 

^ r, ax , . axix ■ x) 
x= - 1+^ c^ — + 2(1 + 7) 

2 _ a^x _ a^xix-xY 



+ 2c'a (/3-e + 27(l + 7)) ^ + 2ae 

X'^ x^ 



+ 2 « (2(1 - /3) + e - 2 7^) ^ ^ " + • (39) 



x"^ 

Setting a = in the solution of fl39|) one can formally get a second-order solution for 
the post-Newtonian equations of light propagation. The merit of this parameter will be 
clear below. 

5.2. Initial value problem 

Let us now solve analytically an initial value problem for the derived equations. For 
initial conditions ( 123|) using the same approach as in P |3] , one gets: 

-Xm = (40) 

c 

Xn = Xo + c{t-to)(T, (41) 

-ipN = (T + mAi{xN), (42) 
c 

XpN = Xn + m {Bi{xn) - Bi{xo)), (43) 

-ippN = cr + mAi(a:;pN) +m^ A2(a;iv), (44) 

,2 



ajppN =xn + m {Bi{xp^) - Bi{xo)) + m {B2{xn) - ^2(2^0)) , (45) 



where 



, . a X Ix X cr) cr \ 
A,{x) = -{1 + ^) — ^ V + - ' (^6) 

\X{X — (7 ■ X) X J 

^ / s / , / cr X (x X cr) , , A , , 

Six =-1 + 7 ^ '- + alog{x + a-x)\, 47 

\ X — cr ■ X J 

. , , 1 cr ■ X , ^r, a X (x X a) , ^r, a x (x x cr) 
A2(a;) = --ae^a; + 2(l + 7)2— L + (i + ^)2 ^ ^ 



2 ^ x'^ x'^ {x - a- ■ x) x{x- a- ■ xf 

2 



-(1+7)'^ ^+ 2(1 -a + 7) (1 + 7) + a/3- -ae - 

X [X — CT ■ X) \ I J X^ 

1 CT X ( iZJ X O" ) 

^° (1 + 7 - a 7) (1 + 7) - 4 a /3 + 3 a e) (cr ■ a;) 



4 |cr X a:;|2 

i (8(i+^-«^)(i + ^)-4«/3 + 3«e) ^'^^ {jx-b{cT,x)) , (4^ 

4 cr X a; r 
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cr , a X (x X cr) 1 x 

2 



B^ix) = -(1 + 7)' + (1 + 7)'- ^ ^ + 7"e 

x-cr-a; {x - (t ■ x) 4 x 

-\a (8(l + 7)-4/3 + 3e) {1-5{ct,x 

_i (8(i+^-«^)(i + ^)-4a/3 + 3ae) {a ■ x) ^ ""^f {-n - 5{ct,x)) , (49) 
4 cr X xr 



or, alternatively, for Si and B2 

, f cr X (x X cr) , , ,\ , , 

Sia; =- 1 + 7 ^ '--a\og{x-a-x)\, (50) 

Y X — <T ■ as y 

a; =+ 1 + 7 ' +1 + 7 '- ^ -^ + -ae — 

x-a x {x-cr -x) 4 x^ 

-i a (8 (1 + 7) - 4 /3 + 3 e) - 6{ct, x 

4 |cr X a^l \2 

(8(l+7-«7)(l + 7)-4a/3 + 3ae) (o- • x) ""i^^^^p"^^ {tx - 5{ct,x)) . (51) 

With these definitions the solution of (139|) reads 
a;(t) =a;ppN(t) + C(c-^), 

^i(t) = ^ippN(t) + 0(c-6). (52) 

It is easy to check that the solution for coordinate velocity of light ippN satisfies the 
integral ( l36l) . In order to demonstrate this fact, it is important to understand that 
position X in ( !36|) lies on the trajectory of the photon and must be therefore considered 
as aSpN in the post-Newtonian terms and as a:;N in the post-post-Newtonian terms of 

5.3. Vector n in the initial problem 
Using (l44l) one gets 

n = (7 + mC i{xpt^) + Ciix^) + C(c"^), (53) 

where 

^ / N . / N / . / NN / , cr X (a? X cr) 

Ci(a^) = Ai(a;) - cr(o--Ai(a^)) = -(1 + 7)- ^ ^ 



x(x — cr ■ x] 



5 



C2ix) = A2{x) - Ai{x) {cr ■ Ai{x)) ~ ^cr {Ai{x) ■ Ai{x)) - a {cr ■ A2{x)) 

+ ^(7 {(T-A,{X)Y 

1 cr ■ X , . , .r, cr X (x X cr) 

= — ae — — cr X (a? X cr) + (1 + 7)^ 



2 X* x^ {x — cr ■ x) 

^r, cr X (x X cr) 1 , ^r, cr x + cr ■ x 

+ l + 7f 4 - o 1 + 7 

— cr-a:;)^ 2 x^ x — cr ■ x 

1 / / ^ , N < ^ N / , cr X fa? X cr) 

- - (8 (1 + 7 - a 7) (1 + 7) - 4 a /3 + 3 a e) (cr ■ a;) 



4 ' ' ' ' " y ■ I ^ ' ■ ' \ ' Icr X a;l2 
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-\ (8(l + 7-«7)(l + 7)-4a/g + 3ae) "^^J^^^^^ (vr - 5((t, a;)) . (54) 

5.4- Impact parameters 

As we have seen in Sections 14.11 and 15.21 the usual analytical solutions are expressed 
through one of the two following impact parameters: 

= (T X [xq X cr) , (55) 
d = k X {xq X k) = k X {xi x k) , (56) 

where Xq is the initial point in both Cauchy and boundary problems given by fl23l) 
and 0151) . respectively, while Xi is the final position in the boundary problem. Both 
these impact parameters naturally arise in practical calculations of light propagation 
when positions of source and observer are given in some reference system (e.g., in the 
BCRS [lOj). However, these parameters are clearly coordinate-dependent and have 
no profound physical meaning. One can expect that formulas involving these impact 
parameters contain some spurious, non-physical terms obscuring the physical meaning 
of the formulas. As we will see below it is indeed the case. Now, we introduce another 
impact parameter 

d' = lim -xit) X (xit) X - xit) \ = lim cr x ixit) x cr) . (57) 
For a similar impact parameter defined at t — > +oo 

d" = lim -x{t) X (x{t) X - x{t)] = lim i/ x {x{t) x u) , (58) 



t->+oo C \ C 

where u = lim -xit), one has Id'l = \d"\. It is also clear that the angle between 

d' and d" is equal to the full light deflection (see below). Since both d' and d" 
reside at time- like infinity and since the metric under study is asymptotically fiat, these 
parameters can be called coordinate- independent. 

One can show that d' = d" coincides with the impact parameter D introduced, 
e.g., by Eq. (215) of Section 20 of [4j in terms of full energy and angular momentum 
of the photon (see also [Ij for a useful discussion). Indeed, in polar coordinates i^x.ip) 
the Chandrasekhar's impact parameter D = f{x)x'^!f, where lim f{x) = 1. Clearly, 
x"^ (f = \x{t) X x{t)\ and it is obvious that d' = d" = D. Interestingly, this discussion 
allows one to find an exact integral of the equations of motion for a photon in the 
Schwarzschild field. The equations of light propagation f|T3|) in the Schwarzschild metric 
([2]) in harmonic coordinates have an integral 
(I + I 

D = ^ - xi^t) X xit) = const, (59) 

1 — a c 

while for the parametrized post-post-Newtonian equations of motion given by fl39|) one 
has 

D = exp(2{l + -f) a + a (2 (1 - /3) + e - 27^) ^x{t) x x{t) 
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= (l + 2(1 + 7) a + (2(1 + jf + a(2{l-f3) + e- a^) ^ x{t) x x{t) + 0{c~') 

= const. (60) 

The first line of (MJh represents an exact integral of the (approximate) equations of 
motion ( 139|) . In both cases the Chandrasekhar's D is the absolute value of D as given 
above. 

Let us stress that the impact parameter d' is not convenient for practical 
calculations, but we will use it below to understand the physical origin of various terms 
in the formulas describing the light propagation. Therefore, we need to have a relation 
between impact parameters (155|1 . (156|1 . and (1571) . Relation between d' and dg. can be 
derived using the post-Newtonian solution for light propagation given above: 

(in \ 
l + + ^ {xo + (T-Xo)j +0{c~^). (61) 

Relation of d ' and d can be derived using formulas of Section 14.11 

d =d{l + il + ^) — 



(PR 2R J 

- (1 + 7) m fc ^i-^o + R ^ ^^^-4)_ ^g2) 
R 

Now we are ready to proceed to the analysis of the post-post-Newtonian equations of 
light propagation. 

5.5. Total light deflection 

In order to derive the total light deflection, we have to consider the limits of the 
coordinate light velocity i for t — )■ ±00. Using formulas of Section [5^2] one gets 

lim - x(t) = cr, (63) 
lim - x(t) = u 

t->-+oo c 

a X [xn X cr) , 9 cr 

= cr - 2 1 + 7 m , ^ V„ ^ -21 + 7 ^ 

|a3o X cr|^ |a?o X crl 

1 / / N / N r. N 9 f X (iCn X cr) 

- -TT f8(l + 7-a7)(l + 7) - 4a/3 + 3ae) ' 



4 \xi^ X cr|3 

X (a^o X 
|a;o X cr| 



n9 9 / , cr X [xn X cr) , , , 

+ 2(1 + 7)2 {xo + (T- xo) ^ + 0(c-6). (64) 



Therefore, the total light deflection reads 
|cr X i^l = 2 (1 + 7) m ^ -2(1 + 7)^ (xq + cr ■ a^o) 



|a;o X cr| | a;o x cr 



3 



+ ^(8(l + 7-a7)(l+7)- 4af3 + 3ae) vrm^- + (9(0-6) . ^q^^ 

4 |a;o X cr|2 

Eq. ( 165|) defines the sine of the angle of the total light deflection in post-post-Newtonian 
approximation. The first term in (16 5 p is the post-Newtonian expression of total light 
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deflection. The other two terms are the post-post-Newtonian corrections. Using d' 
defined by fl57|) and related to d„ by fl6T]) one can rewrite fl65|) as 

|o- X = 2(1 + 7)- + - (8(l + 7-a7)(l +7) - 4a/3 + 3ae) vr — + C(c-^) . (66) 

This result with a = 1 coincides with Eq. (4) of [B] and also agrees with the results 
of [m [5], |2l [18] in the corresponding limits. It is now clear that the second term in 
the right-hand side of (|65|) "corrects" the main post-Newtonian term converting it to 
2(1 +'j)m/d'. Note that the total light deflection |cr x i/| is a coordinate-independent 
quantity and (l66l) expresses it through coordinate-independent quantities while (l65l) 
does not. 



6. Post-post-Newtonian solution of the boundary problem 

For practical modeling of observations it is not sufficient to consider the initial value 
problem for light propagation. Two-point boundary value problem given by (fT5|) is 
important here. This Section is devoted to a derivation of the post-post-Newtonian 
solution of this boundary problem for (|39|) . 

6.1. Formal expressions 

An iterative solution of (HU]) - (H5]) for the propagation time r = ti — to and unit direction 
a reads: 

CT = R-mk- [Bi{xi) - Bi{xo)] - k ■ [^2(^1) - Biixo)] 

+ |fc X {B.M - B,{xo))f + 0{c-'), (67) 

(7 = k + m^{kx[kx {Bi{xi) - Bi{xo))]) 

+ m^-{kx[kx {B2ixi) - B2ixo))]) 
ri 

+ — (Si(xi) - B^{xq)) x[kx {Bi{xi) - Bi{xo))] 

-^m^^k \kx{B^{x,)-B,{xo))f + 0{c-'). (68) 

These expressions are still implicit since in order to achieve the post-post-Newtonian 
accuracy one should use the post-Newtonian relation between a and A; to represent a in 
Bi appearing in the post-Newtonian terms. That relation can be again obtained from 
(1681) by neglecting all terms of order 0{c~^). On the contrary, in the terms of the order 
of 0{c~'^) in fl67|) and fl68|l one can use the Newtonian relation a = k. 
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6.2. The propagation time ct 

Substituting ( l47l) and ( l49l) into ( !67l) one can derive an explicit formula for the time of 
light propagation as function of the given boundary conditions Xq and Xi: 

CT = R 

xi + xo + R 



pN 
ApN 

ppN 
ppN 



+ (1 + 7) m lo. 




i a (8(1 + 7) - 4/3 + 3e) 6 {xi, xq) 

4 |a;i X ajol 



+ 0(c-^) . (69) 

Here we have used that 6{k,Xo) — 5{k,Xi) = 6{xi,xq). Here and below we classify 
the character of the individual terms by labels N (Newtonian), pN (post-Newtonian), 
ppN (post-post-Newtonian) and ApN (terms that are formally of post-post-Newtonian 
order (9(c~'^), but may numerically become significantly larger than other post-post- 
Newtonian terms; see below). Using \xi x a^ol = Rd where d is the impact parameter 
defined by ( l56i) . and assuming general-relativistic values of all parameters a = /3 = 7 = 
e = 1 one gets the following estimates of the sums of the terms labelled by "ApN" and 
"ppN", respectively (the proofs are given in pTl 122] ): 

|c5r^pN|<2^i?-i^l^ <2^i?, (70) 
d^ [Xi + Xq)^ d^ 

I c- , 15 , 

Ic^TppNl < — TT — . (71) 

These estimates and all estimates we give below are reachable for some values of 
parameters and, in this sense, cannot be improved. From these estimates we can 
conclude that among the post-post-Newtonian terms C(5rApN can become significantly 
larger compared to the other post-post-Newtonian terms. For this reason we will call 
such terms "enhanced" post-post-Newtonian terms. The physical origin and properties 
of the "enhanced" post-post-Newtonian terms will be discussed in Section [7l 

The effect of |c5rppN| for the Sun is less than 3.7 cm for arbitrary boundary 
conditions and can be neglected for any current and planned observations. Therefore, 
the formula for the time of light propagation between two given points can be simplified 
by taking only the relevant terms: 

D , n , ^ 1 xi+xo + R 
CT = R + (l-|-7jm log 

- hl + jfm'- ^ (i?^_(a;,-xo)^)+0(::^)+0(m^). (72) 
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This expression can be written in an elegant form 

a^i + + -R + (1 + 7) m 



cr 



R 



7) m lo^ 



R+{l + j)m \ d J ^ ^ 



(73) 



Xi + Xq 

that has been already derived in ^3] in an inconsistent way (see Section 8.3.1.1 and 
Eq. (8-54) of [13]). As a criterion if the additional post-post-Newtonian term is required 
for a given situation, one can use ( ffOl) giving the upper boundary of the additional term. 



6. 3. Transformation from k to cr 
Substituting (gT]) and (gQ]) into ([MD one gets: 
k 

X\ 



pN 



ApN 



scaline 



ppN 



ppN 



ppN 



(1 7) m ■ 



\Xx X ajg 



+ 



2 

(1 + 7)' 



m k X (xq X Xi) 
{xi -Xo + Rf 



(Xi + Xq) {Xi - Xq - R){xi - Xq + Rf 



\Xi X XqI 



■ m 



+ k X [xq X Xi] 



\Xi X Xq\ 



1 



k 

1 



at — - 

4 i?2 



Xi 



■^0 , 



+ - (8(l-h7-a7)(l-F7) -4a/3-h3ae) 



1 



\Xi X Xq\ 



2R^ (tt - 5{k, x)) + (x? -xl- R^) 5{xi, Xq 



+ 0{c 



(74) 



This formula allows one to compute a for given boundary conditions Xq and Xi. Let us 
estimate the magnitude of the individual terms in f l74p in the angle (5(cr, k) between cr 
and k. This angle can be computed from vector product p = k x cr, and, therefore, the 
term in ( 17^ proportional to k and labelled as "scaling" plays no role. Here and below 
terms proportional to k do not influence the directions in the considered approximation. 
These terms are only necessary to keep the involved vectors to have unit length. Now, we 
represent the vector product p as the sum of three kinds of terms: p = ppN+PApN + PppN 
where each term is the vector product of k and the sum of the correspondingly labelled 
terms in flTl]) . Using 

\k X [k X {xq X Xl)] I = I X {xq X Xl) \ = Rd, (75) 

and general-relativistic values of the parameters a = (3 = 'j = e = l one gets (the proofs 
can be found in [211 [22]): 



IPpnI < -r 



Am 
J 



Xl 
Xl + Xq 



Xq < Xl, 



,Xo> Xl 



< 



Am 



(76) 
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IPapnI < 16 



2 

m 



4 . X 1 

^?.Tn , (77) 



{Xi + Xq) 



2 , a;o < f Xi or xq > Xi 



, , 15 , , 

|PppnI<^^^- (78) 

Note that p^^ and Papn themselves as well as their estimates are not continuous for 
Xi ^ Xq since in this limit an infinitely small change of a^i leads to big changes in k. 
Discontinuity of the same origin appears for many other terms. The limit Xi — )■ Xq and 
the corresponding discontinuity have, clearly, no physical importance. 

We see that among terms of order 0{m?) only IpapnI cannot be estimated as 
const X im? /d?. The sum of the three other terms labelled as "ppN" can be estimated as 
given by (ITS]) . In most cases these terms can be neglected at the level of 1 /ias. Indeed, it 
is easy to see that |PppN| can exceed 1 yuas only for observations within about 3.3 angular 
radii from the Sun. Accordingly, we obtain a simplified formula for the transformation 
from to cr keeping only the post-Newtonian and "enhanced" post-post-Newtonian 
terms labelled as "pN" and "ApN" in 

cT = k+dS (l-s]-{xi+xo)(l + +o{'^\+ 0(m=^) , (79) 



2 ' ' V RJJ \(P 

Eq. (177|) can be used as a criterion if the post-post-Newtonian term in ( 17^ is necessary 
for a given accuracy and configuration. 

6.4- Transformation from cr to n 

The transformation between n and cr is given by fl53|l - fl5^ . We need, however, to 
express the relativistic terms in fl53|) as functions of k. To this end we note that 
ajpN = Xi + 0{c~^) and xi = Xq + Rk, use f l71|) for cr in Ci^x^^), and get 

n = cr 



pN 

scaling 

ApN 

ppN 

ppN 

ppN 



\ 1 / \ R ( ^ k- Xi\ 

- (1 + 7)mA; x [xq x xi)- 1 H 

I 331 X a^o r V ^1 / 

+ ^ (1 + 7)2 — ( 1 + (3xi -xo- R) (xi -xo + R) 

4 I a;i X a^o r xi \ xi J 

2 Xi -|- Xo , k ■ xi\ R {R"^ — {xi — xo)"^^ 



+ rn^ k X (xq x x 



(i + 7)^r^^^-?^ 1 



Xi X Xq l"^ \ Xi / 2 \ Xi X Xq P 



/i \2 /i k ■ Xi\ 1 

+ (1 + 7)2- 1 + — 

I £Cl X Xq I \ Xi / Xi 

+ 1(1+7)^^^ + (l-^^l±^) (R^-i 

2^ \xiXXq\^ \ xi J \ R J ^ ^ 

1 k ■ Xi 



Xi - Xo)2 
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ppN 



ppN 



\ (8(l + 7-a7)(l + 7) -4a/3 +3ae)^^^ 

4 xf \xix xo\^ 



4 



- ] (8(1 + 7 - a7)(l + 7) - 4a/3+3ae) 



Xi X Xq 



[n - 6{k,xi)) 



(81) 



This expression allows one to compute the difference between the vectors n and cr 
starting from the boundary conditions Xq and Xi. Let us estimate the magnitude of 
the individual terms in ( IHTjl in the angle S{cr, n) between n and cr. This angle can be 
computed from vector product cp = cr x n. Again the term in (IHT]) proportional to k 
and labelled as "scaling" plays no role since cr x k = 0{c~'^). In order to estimate the 
effects of the other terms in f lSTl) . we split (p = (p^^ + (fi^pN + '^ppN similarly as we did 



with p above, take into account that | cr x (fc x (xq x xi)) 
again a = /3 = 7 = e = l and get [211 



Rd + 0{c ), assume 



I'/^pnI = 
Iv'apnI 



2 m 



cr X [k X [xq X xi)] 



R 



Xi X Xq 



1 



k ■ Xi 

Xi 



m 



^2) 



4m^ 



, m 
< 4 — 



cr x\k X (xq X xi)] 

^ 4Xi Xq 



k ■ Xi\ R {Xi+ Xq) R^ — (Xi — XqY 



X-i 



Xi X X, 







Iv'ppnI<^-- 



d? d (Xi + Xq) 

15 rn? 



(83) 

(84) 



Eq. f lM|) shows that the "ppN" terms can attain 1 /xas only if one observes within 
approximately 3.3 angular radii from the Sun. In many cases these terms can be 
neglected. Accordingly, we obtain a simplified formula for the transformation from 
cr to n keeping only the post-Newtonian and "enhanced" post-post-Newtonian terms 
labelled as "pN" and "ApN" in (EI]): 



n 



(T + dT 1 + Txi 



m 



-(1+7)^ 1 + 



R + Xq — Xi 
R + Xq + Xi 

k ■ Xi\ 

Xi 



m 



0(-l+0(^3) 



^5) 



(86) 



Eq. fl83|) can be used as a criterion if the additional post-post-Newtonian term in fl85|) 
is necessary for a given accuracy and configuration. 



6. 5. Transformation from k to n 

Finally, a direct relation between vectors k and n should be derived. To this end, we 
combine ( FMll and ( I8TII to get 



n = k 



pN 



- (1 + 7) m 



k X {xq X Xi) 

Xi (Xi Xo + £Ci ■ Xq) 
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ApN 



scaling 



ppN 



PPN 



ppN 



+ (1 + 7)^m^ 



Xi Xo + Xi- Xq) 
2 



R 



Xi 
2^2 



+ k X (xq X 



1 1 / 1 



2 ajf \xi X a^ol 



- 2 



k ■ x-i 



ppN 



1 

4 
1 



l + 7-a7)(l + 7) -Aa(3 + 3ae 



R 



k ■ xi 



xf 



Xi X X, 







X 







6{xi,xo) 



+ - (8(l + 7-a7)(l + 7) -4a/3 + 3ae 

8 \Xl Xq]" 

+ 0(c"^) . (87) 

This formula allows one to compute the unit coordinate direction of light propagation 
n at the point of reception starting from the positions of the source Xo and the observer 

Xi. 

As in other cases our goal now is to estimate the effect of the individual terms in (187|) 
on the angle 5{k, n) between k and n. This angle can be computed from vector product 
u: = kxn. The term in fl87|l proportional to k and labelled by "scaling" obviously plays 
no role here and can be ignored. For the other terms in a; = cUpN + ^ApN + <^ppN taking 
into account (!75|) and considering the general-relativistic values a = /3 = 7 = e = l one 
gets 



pN| 



1 |fc X (a;o X xi] 
2m — 



m 

< 4 — 



Xq 



I^^^ApnI = 4m 



Xi Xi Xq + Xi- Xq 
2X1+ Xq \k X {Xq X Xi] 



, m 
< 4 — 

d Xi + Xq d 



Xi 



iXi Xq + Xi- Xq) 



m? RxtX 
< 16 — 



# [xi + Xq) 



m 

< 16 — 



xixl 



d^ {xi + xq)"^ 



x^ 

< 16 

- d^ d 



(89) 



or, alternatively, 



QAm"^ R 
|c^Apn|<^^^. 

We give four possible estimates of jcjApN 
situations. Note that the last estimate in 



(90) 

These estimates can be useful in different 
and the estimate in flUUl) cannot be related 



to each other and reflect different properties of IcjapnI as function of multiple variables. 
The effect of all the "ppN" terms in (IHTj) can be estimated as 



\(JL> 



PpN 



15 



(91) 



Again these terms can attain 1 fias only for observations within about 3.3 angular 
radii from the Sun and can be neglected. Accordingly, we obtain a simplified formula 
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for the transformation from k to n keeping only the post-Newtonian and "enhanced" 
post-post-Newtonian terms labelled as "pN" and "ApN" in flHTj) : 

n = k + dP ^l + Px,^^^^+o(^^^+0{m'), (92) 

Let us also note that the post-post-Newtonian term in fl92l) is maximal for sources at 
infinity: 

I I / r I I r < \2 + |fc x (a^o x Xi)| 

I'^^ApNl < hm IcjapnI = hm (1 + 7) m 



= (l + 7)^(l-cos$)2^^, (94) 

where $ = 6{xo,Xi) is the angle between vectors Xq and Xi. Several useful estimates 
of this term are given by (189|) - (!90|) . These estimates can be used as a criterion which 
allows one to decide if the post-post-Newtonian correction is important for a particular 
situation. 



6. 6. Transformation from k to n for stars and quasars 

In principle, the formulas for the boundary problem given above are valid also for stars 
and quasars. However, for sufficiently large xq the formulas could be simplified. It is 
the purpose of this Section to derive the formulas for this case. 

6.6.1. Transformation from k to cr. First, let us show that for stars and quasars the 
approximation 

(T = k (95) 

is valid for an accuracy of 1 /ias. Using estimates f l76|) and fl77|) for the two terms in 

fl79|l one can see that for xq ^ Xi the angle 6{cr, k) can be estimated as 

, , m Xi / m Xt Xo \ , ^ 

6{cT,k)<A-^^ l+4--i^^ . (96) 

a xi + xo \ a a xi + xo / 

Clearly, 6{(t, k) goes to zero for xq — )■ 00. Numerical values of this upper estimate are 
given in Table [2] for xq equal to 1, 10 and 100 pc. Angle S{<t, k) is smaller for stars at 
larger distances. However, for objects with xq < 1 pc the difference between cr and k 
must be explicitly taken into account. From the point of view of the relativistic model 
these objects should be treated in the same way as solar system objects. 

6.6.2. Transformation from cr to n. As soon as we accept the equality of cr and k for 
stars the only relevant step is the transformation between cr and n. This transformation 
in the post-post-Newtonian approximation is given by f lS^ -f lM|) . In the framework of 
the relativistic light defiection model, the distances to stars and quasars are assumed to 
be unknown and so large that they can be considered infinitely large. For such sources 
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Table 2. Numerical values of estimate (|96p in /las for the angle between cr and k due 
to the solar system bodies for various values of xq. 



Xq [pc] 


Sun 


Sun at 45° 


Jupiter 


Saturn 


Uranus 


Neptune 


1 


8.506 


0.056 


0.473 


0.309 


0.212 


0.382 


10 


0.851 


0.006 


0.047 


0.031 


0.021 


0.038 


100 


0.085 


0.001 


0.004 


0.003 


0.002 


0.004 



it is natural to use the observer's position Xi as initial position denoted in ( ]23l) as Xq. 
Therefore, in fl55|) and fl6T|) one should formally replace Xq by Xi. E.g., the impact 
parameter is defined as 

d„ = (T X {xi X cr). (97) 

We can rewrite fl53l)-(l54l) as 



N 
pN 

ApN 

scaling 

ppN 

ppN 

ppN 

ppN 



n 



l + j)m— [1 + 



dl \ xi 
+ (l + 7)'m2— - 1 



dl d„ V xi 
1 9/, n9 cr cr ■ a3i\ 2 

-5™=(i + 7)=5|(n-— ) 

1 2 '^■^Ij 

- -m ae — d^ 

2 Xi 

+ (1 + 7) m — — 1 + 



1 /o /I , \ /I , \ A n . n \ 2 ^cr cr ■ a?! 

I . . nd 



- - (8(l + 7-a7) (1 + 7) -4a/3 + 3ae) m^-f 



J + 7 - 0:7) (1 + 7) - 4a/3 + 3Q;e) -"l^ (tt - (5(cr, x^) 
4 

+ 0(m3), (9J 



where d^ = \da\ = \cr x xi\. Now we need to estimate the effect of the individual 
terms in (1981) on the angle (5(cr, n) between cr and n. This angle can be computed 
from vector product if) = a x n. The term in (1981) proportional to cr and labelled 
as "scaling" obviously plays no role and can be ignored. For the other terms in 
V' = V'pN + V'ApN + V'ppN taking into account that \(t x d„\ = d„ and considering 
the general-relativistic values a = /3 = 7 = e = lwe get [211 ES] 

|cr X dg-l ( ^ <7 ■ x\\ ^ ^ m 

Xi 



IV'pnI = 2"^(1 + — j<4-, (99) 



IV.^„Nl=W^Hi(l + ^)^16^^i, (100) 

d% da\ X\ J d% da 
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Table 3. Numerical values of the analytical upper estimates of the post-post- 
Newtonian terms of order of 0{m? /d) in and of order 0{m? /dP) in ([71)) . ([5T|) . 
dSZl), and dMl)- 





Sun 


Sun at 45° 


Jupiter 


Saturn 


Uranus 


Neptune 


Ic^TppNl [10-6 m] 


36906.0 


242.9 


0.328 


0.036 


0.002 


0.003 


PppnI, IV'ppnI, ] 
I'^ppnI, IV'ppnI 


10937.4 


0.474 


0.945 


0.120 


0.016 


0.023 



Table 4. Maximal numerical value (|94l) of the "enhanced" post-post-Newtonian term 
in ((92|) for the solar system bodies with parameters given in Table [TJ 



Sun 


Sun at 45° 


Jupiter 


Saturn 


Uranus 


Neptune 


max WApN [A^as] 3192.8 


0.663 X IQ-^ 


16.11 


4.42 


2.58 


5.83 



, , , 15 . s 

IV'ppnI< yvr — . (101) 

The estimate shows that the "ppN" terms can be neglected at the level of 1 /las except 
for the observations within about 3.3 angular radii from the Sun. Omitting these terms 
one gets an expression valid at the level of 1 /las in all other cases: 

Cm?' \ 
— \ + 0(m=^) , (102) 

Q=-il+,)-{l + —l). (103) 

This coincides with ( I92l) -( l93l) and with ( l85l) - (!86l) for xq — oo. This formula together 
with (T = k can be applied for sources at distances larger than 1 pc to attain the 
accuracy of 1 /las. Alternatively, Eqs. ( !92|) -(!93l) can be used for the same purpose giving 
slightly better accuracy for very close stars. However, distance information (parallax) 
is necessary to use (!92|) - (!93|) . 

6.7. Numerical estimates and Monte- Carlo simulations 

Table [3] contains numerical values of the "regular" post-post-Newtonian terms of order 
0{m^/d) in ([69D and of order 0{m^ / d'^) m ([HI]), ([87D, and ([98]). The analytical 

estimates are given by ( 1711) . ( ITHl) . f l8^ . ( 19T|) . and f llOip . respectively. One can see that 
at the level of 10 cm in distances and 1 /las in angles these terms are irrelevant except 
for observations within 3.3 angular radii from the Sun. 

A series of additional Monte-Carlo simulations using randomly chosen boundary 
conditions has been performed to verify the given estimates of the post-post-Newtonian 
terms numerically. The results of these simulations fully confirm all our estimates. 

Using estimate (IMj) and the parameters of the solar system bodies given in Table [1] 
one can compute the maximal values of the "enhanced" post-post-Newtonian term in the 
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transformation from k to n. For grazing rays one can apply cos'l' — —1, while for the 
Sun at 45° one can apply cos $ ~ — 1/ \/2. The results are shown in Table |H Comparing 
these values with those in the last line of Table [H one sees that the "enhanced" post- 
post-Newtonian terms match the error of the standard post-Newtonian formula. The 
deviation for a grazing ray to the Sun is a few /las and originates from the post-post- 
Newtonian terms neglected in fl92|) . 

Vector n computed using fl92l) can be denoted as n'^j^f. The numerical validity of n'^j^ 
can be confirmed by direct comparisons of n'^^ and vector n computed using numerical 
integrations of the geodetic equations as discussed in Section 14.21 For example, the 
results for Jupiter show that the error of n'^^ does not exceed 0.04 /las. The origin of 
this small deviation is well understood and will be discussed elsewhere. 

7. Physical origin of the "enhanced" post-post-Newtonian terms 

We have found above the estimates of various terms in the transformations between 
units vectors cr, n, and k characterizing light propagation. These estimates reveal that 
in each transformation "enhanced" post-post-Newtonian terms exist that can become 
much larger that other "regular" post-post-Newtonian terms. In each case the sum of 
the "regular" post-post-Newtonian terms can be estimated as ^ vr The "enhanced" 
terms can be much larger, being, however, of analytical order m^. In this Section we 
clarify the physical origin of the "enhanced" terms. 

First, let us note that the "enhanced" post-post-Newtonian terms in (^^, ( 15T]) . 
f l87|) . and f l98|) contain only parameter 7. It is clear that these terms come from 
the post-Newtonian terms in the metric and in the equations of motion (parameter 
a does not appear in these terms; see Section IS.l.Sp . Therefore, their origin is the 
formal second-order (post-post-Newtonian) solution of the first-order (post-Newtonian) 
equations given by the first line of (15^ . 

Now let us demonstrate that the "enhanced" terms result from an inadequate choice 
of impact parameters d or in the standard post-Newtonian formulas. Indeed, we can 
demonstrate that the "enhanced" terms disappear if the light deflection formulas are 
expressed through the coordinate- independent impact parameter d' defined by fl57p . 
Eqs. (179])-([80]), (185])-([86D, ([92])-(l93]), and f lT02|) - ffT03|) can be written as 



a =k+d'S' +0 i^—j + 0{m') , (104) 

n =CT + d'T' + o(^^^ +0{m^) , (106) 

r= -(1 + 7)^(1 + —), (107) 

n =k + d'P' + o(^j +0(m^) , (108) 
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(109) 



(110) 



(111) 



respectively. Therefore, in each case the "enhanced" post-post-Newtonian terms only 
correct the post-Newtonian terms that use inadequate impact parameter. Let us stress, 
however, that for practical calculations ([79]) -([80]), ([85]) -([86]), ([92]) -([93]), and ffT02]) -ffT03D 
are more convenient. 

8. Summary and concluding remarks 

In this paper the numerical accuracy of the post-Newtonian and post-post-Newtonian 
formulas for light propagation in the parametrized Schwarzschild field has been 
investigated. Analytical formulas have been compared with high-accuracy numerical 
integrations of the geodetic equations. In this way we demonstrate that the standard 
post-Newtonian formulas for the boundary problem (light propagation between two 
given points) cannot be used at the accuracy level of 1 /ias for observations performed 
by an observer situated within the solar system. The error of the standard formula may 
attain ~16 /ias. Detailed analysis has shown that the error is of post-post-Newtonian 
order 0{m'^). On the other hand, the post-post-Newtonian terms are often thought to 
be of order m? /d? and can be estimated to be much smaller than 1 /xas in this case. 
To clarify this contradiction we have derived and investigated the explicit analytical 
post-post-Newtonian solution for the light propagation. For each individual term in 
the relevant formulas exact analytical upper estimates have been found. It turns out 
that in each case there exist post-post-Newtonian terms that can become much larger 
than the other ones and cannot be estimated as const x m?/(P. We call these terms 
"enhanced" post-post-Newtonian terms. These terms depend only on 7 and come from 
the second-order solution of the post-Newtonian equations of light propagation (Eq. 
(39) with a = 0). For this reason one could argue that the "enhanced" post-post- 
Newtonian terms should not be called "post-post-Newtonian", but better "m^-terms" 
or similarly. The physical origin of the "enhanced" terms is discussed in the previous 
Section. The derived analytical solution shows that no "regular" post-post-Newtonian 
terms are relevant for the accuracy of 1 /ias in the conditions of planned astrometric 
missions (Gaia, SIM, etc.). Most of the "regular" terms come from the post-post- 
Newtonian terms in the metric tensor. It is not the post-Newtonian equation of light 
propagation (Eq. (39) with a = 0) itself, but the standard analytical way to solve this 
equation that is responsible for the numerical error of 16 /ias mentioned above. 

The compact formulas for the light propagation time and for the transformations 
between directions cr, n and k have been derived. The formulas are given by (173]) . ( 179]) - 



([HO]), ([H5])-([H6]), ([92])-([93]), and ffT02l)-ffT03l). These formulas contain only terms (both 
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post-Newtonian and post-post-Newtonian) that are numerically relevant at the level of 
10 cm for the Shapiro delay and 1 /ias for the directions for any observer situated in the 
solar system and not observing closer than 3.3 angular radii of the Sun. 

Let us finally note that the post-post-Newtonian term in (p2l) - (!93l) is closely related 
to the gravitational lens formula. Here we only note that all the formulas for the 
boundary problem given in this paper are not valid for d = {d always appear in the 
denominators of these formulas). On the other hand, the standard post-Newtonian 
lens equation successfully treats this case known as the Einstein ring solution. The 
relation between the lens approximation and the standard post-Newtonian expansion is 
a different topic which will be considered in a subsequent paper. 
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